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Abstract 

Langevin equations are used to model many processes of physical interest, 
including low-energy nuclear collisions. In this paper we develop a general 
method for computing probabilities of very rare events (e.g. small fusion cross- 
sections) for processes described by Langevin dynamics. As we demonstrate 
with numerical examples as well as an exactly solvable model, our method 
can converge to the desired answer at a rate which is orders of magnitude 
faster than that achieved with direct simulations of the process in question. 

PACS: 24.60.Ky 25.70.Jj, 25.70.-z, 82.20.Fd 
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INTRODUCTION 



Langevin methods offer a powerful tool for the numerical study of low-energy nuclear 
processes, such as fission and heavy-ion fusion. The evolution of nuclei during such events 
is typically described using a few collective degrees of freedom, evolving under both conser- 
vative and non-conservative forces. The latter, arising from the coupling of the collective 
variables to the intrinsic nucleonic degrees of freedom, can be modeled by a noisy and a dis- 
sipative term in a Langevin description of the collective motion Once such a stochastic 
equation of motion has been written down, it is straightforward to numerically simulate the 
process in question, using a random number generator to supply the noise. By repeating the 
simulation - with different sequences of random numbers - one obtains independent real- 
izations of the process in question, reflecting the statistical distribution of events occurring 
during an experiment. 

The "direct simulation" method outlined above becomes impractical when studying rare 
outcomes. For instance, if we are interested in computing the very small cross-section for 
the fusion of two heavy nuclei, then the vast majority of realizations will end with the nuclei 
flying apart, and the number of simulations required to obtain even a handful of fusion 
events may well be prohibitively large. The purpose of the present paper is to propose a 
general strategy for getting around this problem. 

The basic idea which we shall present is essentially a dynamical variant of importance 
sampling, which amounts to gaining information about one probability distribution (a "tar- 
get" distribution, T), by choosing randomly from another (a "sampling" distribution, S) 
defined on the same space, and then biasing - assigning weights to - the points sampled. 
The weights are assigned in such a way that the biased average of any quantity, over 
points drawn independently from S, and the unbiased average of that same quantity over 

points drawn from T, converge to the same value in the limit of infinitely many sam- 
ples, oo. If the biased average converges faster with N than the unbiased one, then 
importance sampling becomes a practical tool for increasing the efficiency of the numerical 
estimation of the desired average. 

In our case, we are interested in Langevin trajectories describing (for instance) the col- 
lision of two heavy nuclei, with a very small probability for fusion. Our target ensemble, T, 
is then the statistical distribution of all such trajectories with, say, a given initial center-of- 
mass energy and impact parameter. The probability of fusion which we wish to compute 
is defined with respect to this ensemble of trajectories. Our sampling ensemble, S, is the 
distribution of trajectories evolving - from the same initial conditions - under a modified 
Langevin equation, which by design is far more likely to result in fusion. The scheme which 
we propose involves running a number of simulations with the modified equation of mo- 
tion (thus obtaining fusion events with good statistics), then computing the desired fusion 
probability from this data, by biasing each trajectory. 

The organization of this paper is as follows. In Section |, we derive the method which 
we propose in this paper. This section begins with a general discussion of importance sam- 
pling, before focusing specifically on its application in the context of Langevin trajectories. 
In Section || we briefly touch on a number of practical issues associated with the actual 
implementation of the method. Section ^ presents numerical results, in which we illustrate 
the efficiency of the method using a schematic model of nuclear collisions. In Section |^ we 
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analyze how the method behaves for the case of an exactly solvable model, and we conclude 
in Section 0. 

We should stress at the outset that, while we shall present our result within the specific 
context of computing nuclear fusion probabilities, the method we propose can in principle be 
applied quite generally, whenever one is interested in rare outcomes of processes described 
by Langevin dynamics. For a different approach to studying rare events (developed in the 
context of chemical transitions), see the recent work of Chandler and collaborators 0. 



I. THEORY 



A. Importance Sampling 



Importance sampling is based on a very simple idea, embodied by Eq.^ below. Suppose 
we have some space (^-space) on which are defined two normalized probability distributions, 
Ps{C) andpT(C); corresponding to "sampling" and "target" ensembles, S" and T. Supposing 
furthermore that ps{0 > whenever pt(C) > 0, let us introduce a biasing function 



w{0 



MO 
PsiCY 



(1) 



defined at all points ( for which ps{C) > 0. Now let {0)s and {0)t denote the averages of 
some observable 0{() over the two distributions: 



{Oh 



i = S,T. 



(2) 



If we are interested in computing {0)t, then we can do so by repeatedly sampling from the 
target ensemble T: 



N 



{0)t = ^^Jl/N) E oiC 



(3) 



where (i, Cj, ■ ■ ■ is a sequence of points sampled independently from T. However, it follows 
from Eqs.[l| and ^ above that we can equally well express the desired average as: 



N 



{0)t = {wO)s = jiin (1/iV) E ^(Cf ) OiC 



(4) 



ra=l 



where Cf , Cf ^ ■ ■ ■ is a sequence of points sampled independently from S. Thus, provided 
we can compute w{() and 0{Q for any Eq.H gives us a prescription for determining the 
average of O over the target ensemble T, using points drawn from the sampling ensemble S. 
This prescription becomes a useful tool if a sampling distribution can be chosen for which 
the rate of convergence with the number of samples (A^) is faster when using Eq.^, than 
when sampling directly from T (Eq.^). 

Let us now consider how we might apply importance sampling to the problem which 
interests us, namely, computing the probability of fusion for two colliding nuclei. For sim- 
plicity of presentation, we assume for the moment that the collision process is described by 
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the evolution of a single collective degree of freedom (e.g. a distance coordinate between the 
two nuclei), obeying a Langevin equation of the form 

^=vo{x)+m. (5) 

Here, x is the collective variable, vo{x) is a "drift" term which embodies the deterministic 
forces - both conservative and dissipative - acting on the collective degree of freedom, and 
^(t) is a stochastic, white noise term: 

mi{t + s)) = DS{s), (6) 

where (• • •) denotes an average over realizations of i{t), and > is a diffusion constant. 
[More generally, the number of variables required to specify the instantaneous state of the 
system will be greater than one: x — >• x = (xi, ■ ■ ■ , Xd)- In the case of overdamped motion, 
this vector specifies the instantaneous configuration of the colliding nuclei. For evolution in 
which inertial effects are important, the vector x will include both configurational variables 
and their associated momenta. See Section |T| for elaboration of these and other points.] 

Let us assume initial conditions corresponding to a particular value x^ for the collective 
variable, and suppose we are interested in the evolution of the colliding nuclei over a time 
interval < t < r. Then a single realization of this process is described by a trajectory x{t), 
< t < T, satisfying x(0) = x^, which obeys Eq.|] for a given realization of the noise term 
^(t). Therefore, the statistical ensemble of realizations of ^{t) - together with EqJ^ and the 
initial condition x(0) = x^ - defines a statistical ensemble of trajectories x{t). By simulating 
such a process numerically - using a random number generator to provide the noise term 
- we are, effectively, sampling randomly from this ensemble. Let T ("target") denote the 
ensemble of trajectories thus defined: 

T={x{t)\x = vo + i,x{0)=x'^}, (7) 

with the time interval < t < r left implicit, here and henceforth. 

Given the above Langevin equation (Eq.|^), along with the initial condition we are 
interested in computing the probability of fusion. Here we will take "fusion" to mean simply 
that the final point of the trajectory, x(r), falls within a given region R along the x-axis: 

x(r) E R fusion (8) 

x(r) ^ R no fusion. (9) 

Now introduce a functional B[x(t)] which is equal to 1 if the trajectory yields a fusion event, 
and otherwise.[] Then the probability of fusion, Pfus, is just the average of this functional 
over the ensemble T, defined by Eq.^ above: 

Pfus = (e)T. (10) 



^ Of course, in our case, the functional 0[2;(t)] is really just a function of the final point, x(r). 
More generally, however, we could have introduced a criterion for "fusion" which depends on the 
entire trajectory, in which case 0[2;(t)] would be a genuine functional. 
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We can compute Pfus by sampling randomly from the ensemble T (i.e. repeatedly simu- 
lating trajectories evolving under Eq.^, and counting the number of fusion events: 



1 ^ 



^fus = 4^^ = ^E0[^nW], (11) 



n= 



where x'^{t) is the ra'th of independent simulations of the process, all launched from 
x". However, since the number of simulations needed to compute Pfus to a desired accuracy 
grows as the inverse of the fusion probability itself (see Eg .^91) , this "direct sampling" method 
becomes impractical for very small values of Pfus- Let us therefore consider a modified version 
of Eq.|: 

dx 

— =vo{x) + Av{x)+at), (12) 

where the modification consists of adding an extra drift term, Av, chosen to greatly increase 
the probability of simulating a fusion event. Let S (for "sampling") denote the statistical 
ensemble of trajectories x{t) corresponding to this equation of motion, with the same initial 
conditions as above: 

S = [x{t) \ x = vo + Av + x(0) = (13) 

Now, for a given trajectory x{t) satisfying a;(0) = x^, let PT[x{t)] denote the probability 
density that we will obtain this particular trajectory, under the original Langevin equation 
(Eq.D; and let Ps[x{t)] denote the probability density for obtaining this trajectory under 
the modified equation (Eq.|l^). Finally, let w denote the ratio of these two densities^: 

^[,(t)] ^ ^44^)]. (14) 

L ^ Ps[x{t)] ^ ' 

Then, by Eq.^, we can express the fusion probability Pfus (defined with respect to the original 
Langevin equation) as: 

1 ^ 

Pfus = {wQ)s = lim - E ■ 0[^n(i)], (15) 

n=l 

where a;;^(t) is the n'th of independent realizations obeying the modified Langevin equa- 
tion. This means that, if we know how to compute w for any trajectory x{t), then we can 
estimate Pfus by running N simulations with the modified Langevin equation, then adding 
together the weights w for all those trajectories which lead to fusion, and dividing by the 
total number of simulations, N. In the limit — > oo, this estimate converges to the correct 
value of Pfus- 

In the following section we argue that there indeed exists an explicit, easily computable 
expression for w[x(t)], which renders practical the importance sampling strategy outlined 
above. 



^ Note that while the values oips[x{t)] and PT[x{t)] depend on the measure chosen on path space, 
the ratio w[x{t)] is independent of that measure. 
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B. Statistical distributions of Langevin trajectories 



The original and modified Langevin equations, Eqs.^ and |T^, can be represented by the 
generic equation 

-^ = v{x)+m, (16) 

where v = Vq in one case, and v = Vq + Av in the other. As before, given some initial 
conditions x{0) = let x{t) denote the trajectory evolving from those initial conditions, 
for a particular realization of the noise term. We are interested in the probability density 
p[x(t)] for obtaining a particular trajectory x(t). To introduce a measure on the space of all 
possible trajectories, we discretize the trajectory. Thus, let ■ ■ ■ ,x^) denote the set 

of points specifying the state of the system, after time intervals 6t = r/M: 

x'^ = x{m5t) , m = 0,---,M. (17) 

We can then ask for the probability density p(a;^, ■ ■ ■ , x^^|a;°) that the trajectory passes 
through the sequence of points (x^, x^, ■ ■ ■ , x^^) at times 5t, 25t, ■ ■ ■ , r, given the initial 
point Note that this is a probability distribution in M-dimensional (x^, ■ ■ ■ , x*^)- 

space, with x° acting as a parameter of the distribution. Let us introduce the notation 
X = (x*^, x^, ■ ■ ■ , x*^) = (x°, Y). Then an explicit expression for p is given by 

p(Y|x°) = {21x051)-"^'^ exp -A(X), (18a) 

where 

1 2 
.4(X) = — — [a^™^^ - x"^ - v{x'^)6t\ . (18b) 

Eq. |18a| is strictly speaking valid only in the limit M ^ oo (with r fixed), although it 
constitutes a good approximation if V D5t is small in comparison with the length scale set 
by variations in v{x). 

Eq.|l8] is a general expression for the probability density of obtaining a particular dis- 
cretized trajectory X, launched from x''. Now let pt(Y|x°) and p5(Y|x°) denote this prob- 
ability density, for the specific Langevin processes described by Eqs.| and |12|, respectively. 
Then the ratio between these two probability densities is given by: 

^(x) ^ = ^19^) 

Ps[Y\x^) 

AA = At~ As, (19b) 

where At and As are computed with Eq. |18b| , using v = vq and v = vo + Ai;, respectively. If 
we now explicitly write out the expression for in terms of (x*^, x^, ■ ■ ■ , x^), and consider 
the limit M ^ oo (with r fixed), then we arrive at the following result for our biasing 
function, expressed in terms of an integral over a trajectory x{t) rather than a sum over 
points along the trajectory: 

w[x(t)] = exp-AA (20a) 
AA[x{t)] = ^j^^dt(^-v,- Iat;^ A^. (20b) 

Here, dx/dt, vq and Av are evaluated along the trajectory x(t). 
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C. Computing probabilities of rare events 



Combining Eqs.|l^ and ^ we now have an expression which allows us, in principle, to 
compute the probability for fusion - defined with respect to the original equation of motion 
(Eq.^ - by running independent simulations with the modified equation of motion (Eq.|l^): 

1 ^ 

^^^^ = aS^'oo N ^ ^f'^"^^^] ■ ^""P -^^^^nit)]- (21) 
n=l 

Here, x^{t) is the trajectory generated during the ra'th simulation, using the modified 
Langevin equation; AA is computed for each trajectory^; and 6 is, as before, equal to 
one or zero, depending on whether or not fusion occurred. 

For a given set of original and modified Langevin equations, and for a large number N 
of trajectories simulated under the modified equations, we thus have the following estimator 
for Pfus: 

Pf,s = p!^^ = IT^ E 0" exp -AA^, (22) 

n=l 

where 9^ = 6[xf(t)], and A An = AA[x^(t)]. The estimator P^^^ converges to the true 
value Pfus in the limit of infinitely many trajectories: limjy^oo -Pfu^^ = -Pfus- 

D. Efficiency analysis 

Having derived an estimator for Pfus based on the idea of importance sampling, we now 
consider the question of efficiency. In particular, we establish a specific measure of "how 
much we gain" by using importance sampling, with a given choice of Av{x). 

The validity of Eq.^ does not depend on the form of Av{x). Therefore, for any additional 
drift term Av, there will be some threshold value N^^ such that Pf^^ provides a "good" 
estimate of Pfug for N > N^^. That is, N^^ is the number of trajectories which we need to 
simulate (with the modified Langevin equation), in order to determine Pfus to some desired 
accuracy, using the method outlined above. Of course, A^^^, can depend strongly on the form 
of Av{x). We can thus compare the efficiency of estimating Pfus, for different drift terms. In 
particular - since the special case Av = is equivalent to computing Pfus using the original 
Langevin equation - let us define the efficiency gain, E^^, associated with a given Av{x), 
as follows: 

N* 

Et - |f . (23) 

The numerator is just the number of trajectories needed to accurately estimate Pfus by 
running simulations with the original Langevin equation [Av = 0); the denominator is the 
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While Eq.20b gives a well-defined expression for A^, in practice it is more convenient to compute 



AA using tlie method described in Section ^ below; see Eqs.^ to 34 
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number needed using Eq.|22|, for a given Av{x). Thus0 , is the factor by which we reduce 
the computational effort, by making use of importance samphng - again, for a given Av{x). 

Let us derive an expression for E^^ in terms of quantities extracted directly from nu- 
merical simulations. For a given additional drift term Av, let us define 



f[x{t)] = w[x{t)] ■ Q[x{t)] = e exp -AA. 



(24) 



Our method of computing Pfus then amounts to computing the average of /, by sampling 
trajectories x{t) from the ensemble S: {f)s = -Pfus — Pf^^ = {^/N)J2nfn, where = 
f[x^{t)]. The statistical error in our result - the expected amount by which P^^^ differs 
from Pfus - is given by the usual formula for the standard deviation of the mean: 



a 



(s) 

Ptus 



w 



where a? is the variance of the quantity / over the sampling ensemble. 



{f)s 



(/)!• 



(25) 



(26) 



If we want to compute Pfus to a desired relative accuracy r - in the sense that we want the 
ratio (Tpj^^/Pfus (expected error / desired average) to fall below r - then we get the following 
expression for the minimum number of trajectories needed: 



N: 



a 



Ad 



' fus 



(27) 



In other words, if we simulate more than this many trajectories, then we can expect the 
statistical error in our final estimate of Pfus to be no greater than r times Pfus. 

In the case of direct sampling from the target ensemble T, we estimate Pfus by computing 
the average of G. The expected statistical error in this average is just 



ST) 



N 



1/2 



-Pfus(l — -Pfus) 

N 



1/2 



(28) 



making use of the fact that 0^ = 0. By setting this expected statistical error equal to rPfus, 
we get 



AT* _ '^^^^ 
^^0 - ^2 



1 



P 



fus 



r^Pfus' 



(29) 



for Pfus -C 1. (If we simulate this many trajectories, then the expected number of fusion 
events is and the expected statistical error in this number is 1/r.) 

Combining Eqs.^ and we get the following result for the efficiency gain of our 
importance sampling method, for a particular choice of Av{x): 



We are assuming that simulating a single trajectory with the modified Langevin equation takes 
as much time as simulating one with the original Langevin equation. 
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E 



fusi 



Av 



Av 



-r fus 



{f)s 
{P)s-{f) 



(30) 



Eq.^ gives the efficiency gain of importance sampling, with a particular choice of Aw, 
in terms of averages which can be estimated from simulations performed under the modified 
Langevin equation alone (i.e. sampling only from S, not T). An expression for in terms 
of averages estimated using both the original and modified equations of motion, is: 



E 



G 
Av 



AT) 
(S) 



(e)T-(e)| 

{p)s-{m' 



(31) 



We will use these results in Sections |T| and |^ below, to compute the efficiency gain of 
our importance sampling method for particular examples. 



II. PRACTICAL MATTERS 

In this section, we discuss a number of practical issues related to the actual implemen- 
tation of our method. 



At the end of Section [B, we obtained an expression for AA as a functional of the 



trajectory x{t). Since x{t) satisfies Eq.|T2|, we can rewrite Eq. |20 b| as: 

AA = ^ r dt {2i + Av)Av. (32) 
2D Jo 

This expression for AA lends itself to a convenient implementation of our method, as follows. 
When simulating a given trajectory x{t) evolving under Eq.|l^, we simultaneously integrate 
the following equation of motion for a new variable y{t), satisfying the initial condition 
y{Q) = 0: 

^ = ^(2e + A.), (33) 

for the same realization of the noise term ^(t). (Note that this equation is coupled to the 
equation of motion for x, since Av in general depends on x.) Eq.|3^ then implies that 

AA = yiT). (34) 

Thus, at the end of the simulation, we use x(r) to determine whether or not fusion has 
occurred, and if so, then we take AA = yij) when assigning the bias e~^^ to this event. 

• Often (see for instance Section |TTl| below), the evolution of our system is such that, 
once a trajectory x{t) enters the region R which defines fusion, its chance for subsequently 
escaping that region is negligible: R effectively possesses an absorbing boundary. If this is 
true for both the original and modified evolution, it becomes convenient to define Af to be 
zero everywhere within R. Then, if a trajectory x{t) (evolving under the modified Langevin 
equation) crosses into R at some time r' < r, we can stop the simulation at that point in 
time, and take B = 1, AA = y{T'). This saves time, by eliminating the need to continue 
with the simulation. 
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• We have, to this point, assumed that the stochastic noise ^(t) is independent of x. 
That is, D =const. More generally, we might have a diffusion coefficient which depends on 
the instantaneous configuration of the system: D = D{x), i.e. 

{i{mt + s)),(^t)=, = Dix)6{s). (35) 

In this case Eq.^ becomes 

rM-i -1/2I 
p{Y\x^) = lll[27rD{x"')6t\ exp -A(X) (36) 

lm=0 J 

A{X) = —y — r - - v{x"')6t\ . (37) 

Eq.|T9| remains unchanged (since the factor multiplying in Eq.^ is the same for ps and 
Pt)', and Eq.|2y changes only in that 1/D is brought inside the integral in Eq. |20b| , with D 



evaluated along the trajectory x{t). When implementing the method using the additional 
variable y{t), the only difference is that D is evaluated along x{t) rather than being a 
constant, in Eq.|33|. 

• Let us now drop the assumption that the system evolves in one dimension. The state 
of the system is now described by a vector x = (xi, ■ ■ ■ , xa), evolving under a set of coupled 
Langevin equations, 

dx ' ^ 

-^ = v,{x)+Ut) , t = l,---,d, (38) 

where v = Vq (original equations of motion) ot v = vq + Av (modified equations of mo- 
tion). The diffusion coefficient D becomes a symmetric matrix whose elements reflect the 
correlations between the different components of the stochastic force: 

{mUt + s)) = D,,5is). (39) 

[For simplicity, we assume that this matrix is a constant. The generalization to Dij = Dij{x) 
is as straightforward as in the one-dimensional case.] 

The simplest case of mult i- dimensional evolution occurs when the components C,i are 
mutually uncorrelated. Then D is a diagonal matrix, and the generalization of Eqj3^ is 
given by: 

I- E |^(2|. + A„,). (40) 

The sum is taken over values of i for which Da 7^ 0, corresponding to those directions along 
which there is a non-zero stochastic force. Along those directions for which Da = 0, we 
must have Avi = 0.0 

If D is not diagonal, then Eq.EO generalizes to the following evolution equation for y: 



^ To see this, note that if Da = for a particular value of i, then the equation of motion describing 
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| = l(2f+A.rZ^-^At;. (41) 

Eq.^ implicitly assumes that D is invertible, i.e. det(il') 7^ 0. If this is not the case, then - 
first - we must make sure that the projection of At? onto the subspace spanned by the null 
eigenvectors of D is zero. (See the comments following Eq.^) Assuming this condition is 
satisfied, we can view Eq.^ as pertaining only to the subspace spanned by the non-zero 
eigenvectors of D. 

The results discussed in this generalization from one- dimensional to mult i- dimensional 
evolution are based on the following generalization of Eq.|T8|: 

p(Y|x°) = [(27r5t)'^'det(Z))]"*'^^%xp-A(X) (42a) 

1 M-l 

A{X) = - - v{x^)5tY D-\or+^ -x^- v{x^)5t]. (42b) 

As with Eq.^ above, these equations pertain to the d'{< (i)-dimensional subspace spanned 
by the non-zero eigenvectors of D. 

There is an example of multi-dimensional evolution worthy of particular mention. This 
is the case in which inertial effects are present, i.e. the evolution of the system is not over- 
damped. The evolution then occurs in the phase space of the system, and the equations of 
motion for x = {q,p) are typically of the form: 

f=/-V (43a) 
^=F-TI-'p + C (43b) 

Here, g* is a vector of variables specifying the configuration of the system; p is the vector of 
associated momenta; / is an inertia tensor; F is the vector of conservative forces acting on 
the system; F is a friction tensor; and ^ is the vector of stochastic forces, whose associated 
diffusion tensor D is related to F by a fluctuation-dissipation relation. (Typically, /, F, F, 
and D are all functions of q.) In this case, the equations of motion for q are deterministic 
(Eq. [43a| ), therefore any additional drift terms must appear only in the momentum equations 



(Eq. [43b| ), as an additional force AF. The equation of motion for y{t), Eq.^, then pertains 
only to momentum space (the j>-subspace of phase space). 

• Finally, it is often the situation that the dissipative and stochastic forces acting on the 
collective degrees of freedom depend on the "temperature" of the bi-nuclear system. This is 



the evolution of Xi is deterministic (^j = 0) . If we now modify that equation by adding an additional 
non-zero term Avi, then any trajectory x(t) obeying the modified equations of motion will not be 
a solution of the original equations - for any realization of the noise term ^(t) - and vice- versa 
(unless Avi happens to be zero exactly along the trajectory). This violates the condition stated 
before Eq|l|: in order for the importance sampling to be valid, our modified equations of motion 
must be capable of generating any trajectory which might be generated by the original equations 
of motion. 
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another way of saying that these forces, at time t, depend on the total amount of collective 
energy which has been dissipated up to that time. Since the energy dissipated, as a function 
of time, will differ from one realization to the next, it seems we are faced with "memory- 
dependent" forces, i.e. forces which, at time t, depend not only on the instantaneous state of 
the system, but also on its history, up to time t. An easy way to deal with this situation is 
simply to expand our list of variables x to include a new member, Xd+i = -Ediss, denoting the 
collective energy dissipated. This new variable is initialized at zero, and evolves under an 
evolution equation which depends on the model used to describe the colliding nuclei. (For 
instance, 

^ = (/-Vr[r/-V-a (44) 

in the case of collective evolution described by Eq.^.) 

With the addition of this new variable, i.e. a; — >• {xi, ■ ■ ■ ,Xd, Xd+i), we now again have a 
set of (coupled) Langevin equations, in which the drift and stochastic terms depend only on 
the instantaneous state of the system, x. We can thus apply the method proposed in this 
paper, without further modification. 



III. NUMERICAL RESULTS 

In this Section we describe numerical experiments which we have carried out to test our 
method, using a simphfied model of heavy ion collisions introduced by Swi§tecki [Q. This 
model was previously studied by Aguiar et al in 1990 ||^, using Langevin simulations. For 
our example, we considered the collision of two ^°°Zr nuclei. In this mass-symmetric case - 
for this simple model - the shape of the system is defined by two equal spheres connected 
by a cylinder. There are two macroscopic ("collective") variables parametrizing the shape: 
(1) the relative distance p between the sphere centers, which is the distance s divided by 
the sum of radii of the two spheres: p = s/2R] and (2) the window opening a, which is the 
square of the ratio of the cylinder radius to the radius of the sphere: a = {rcyi/RY- 

After some approximations for the potential, dissipation and kinetic energy terms, one 
obtains the following coupled differential equations for the time evolution of the system (see 
Ref. [| for details): 

d^cr r,d(T - , , 

du 2z/ + 3z/2-a ; 

1 A f ^ 2^ = ^2- 45b 

ar Apia p-^) 



(While Eq. [45a| has been written here as a second-order stochastic differential equation, in 



practice we convert it to two first-order equations - one deterministic, one stochastic - 
by introducing the variable p„ = pda/dr.) Here, the collective coordinates p and a are 
represented by the variables z/ = y/a and a = p"^ — 1; the constant yU is a reduced mass, r 
is a reduced time, X is a constant conservative force, and and ^2 are stochastic forces 
(Gaussian white noise), related to the dissipative terms by a fluctuation-dissipation relation. 
(For the ^°°Zr+^°°Zr collision, we have p = 0.176 and X = 0.677.) The evolution of the 
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colliding nuclei is then represented by a Langevin trajectory in (cr, i/)-space. Fig.|l| depicts 30 
such trajectories, all starting from a configuration of two touching spheres {a = 0,i> = 0), 
with a center-of-mass energy equal to 0.8 MeV above the interaction barrier. This energy 
is about 2.5 MeV below the "extra push" energy, so most of the trajectories (28 of them) 
lead to reseparation of the system (fission), and only two trajectories lead to a compound 
nucleus (fusion). 

From Fig.jl] we have the following picture of the physical process occuring, in the context 
of this simplified model: first the window opening between the two nuclei grows rapidly; then 
around a saddle point, at (ex, u) ~ (0.0, 0.6), the combination of deterministic and stochastic 
forces determines the ultimate fate of the nuclei, either fusion or reseparation; and finally 
the system evolves toward its destiny, with a decreasing in the case of fusion, or increasing 
with reseparation. This suggests that, if we are to add an extra drift term to increase the 
likelihood of fusion, then it would be best to localize such a term in the vicinity of the saddle 
point. We will think of such a term as a force, pushing the system toward fusion. We have 
chosen an additional force along the (negative) a direction, whose strength is a Gaussian 
function of (a, z/), with a peak at (0.0,0.6). This leads to the following modified Langevin 
equations of motion: 



with A an adjustable parameter. Fig.^ schematically shows the region around the saddle 
point, where the additional force pushing the system toward fusion is localized. Two de- 
terministic trajectories (evolving under the original equations of motion, but without the 
stochastic terms) are also shown, to guide the eye. One of these was launched with an energy 
of 1 MeV above the barrier (leading to reseparation), the other one at 5 MeV above the 
barrier (leading to fusion). 

To compare our importance sampling method to direct simulation of the original process, 
we first chose to compute the probability of fusion for trajectories starting with an energy 
0.2 MeV above the barrier. This probability is on the order of 10~^, considerably less than 
that for the case shown in Fig.|l| (0.8 MeV above the barrier). We ran 10^ independent 
trajectories under both the original and the modified Langevin equations, Eqs.^ and 
respectively, and kept a running tally of the probablity of fusion, as computed by the two 
methods (Eq.|ll]for the case of direct simulation, Eq.|2^ for importance sampling). We took 
the strength of the additional force to be A = 0.3/i in these simulations. Fig.^ illustrates the 
difference between the rates of convergence of the two methods. The plots show the estimates 
Pf^^ as computed by each method, as a function of number of trajectories simulated, A^. 
It is clear that the importance sampling (broken line) converges much faster than direct 
simulation (solid line): after about 5000 trajectories, the former has converged very close to 
its asymptotic value, whereas the latter is still "jumping around" after 10000 trajectories. 

The sawtooth pattern exhibited by the direct simulation estimate is typical of the sit- 
uation in which rare events contribute disproportionally to an ensemble average: in this 
particular case, only 15 of the 10000 "direct" trajectories (simulated under Eq.^Sj) lead to 
fusion; each of these events causes a sudden jump in the Pf^\ followed by a gradual de- 




(46a) 



du 2v + 3z/2 - o- 



(46b) 



dr 4:u{a + z/^) 
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cline as non-fusion events accumulate. By contrast, under the modified Langevin equations, 



Eq.|46|, about 25% of the 10000 trajectories lead to fusion, resulting in smoother and faster 
convergence. 

In Figs.§ and § we show excitation functions - fusion probability plotted against center- 
of-mass energy above the barrier - as computed by both direct simulation and importance 
sampling, with the same additional force as used in Fig.Q. Each point was obtained using 
1000 trajectories, and the result is displayed with error bars, as estimated from the numer- 
ical data. The solid line represents an analytical formula which closely approximates the 
fusion probability over the region shown.[| Again we see that, for approximately the same 
computational effort, our importance sampling method gives significantly better results than 
direct simulation. For the point corresponding to 0.5 MeV above the barrier, the error bar 
in Fig.^is about 5.5 times bigger than that in Fig.^ The efficiency gain of the importance 
sampling approach is therefore about 30 (~ 5.5^, see Eq.^) in this case: we would need to 
launch about 30 x 10'^ trajectories evolving under the original Langevin equation to get the 
same degree of accuracy obtained in Fig.|^ with 10^ trajectories. 

The gain in efficiency becomes more dramatic when we go to very small probabilities. 
To show this we considered the reaction ^^°Pd-|-^^°Pd [fi = 0.174, X = 0.794), for which 
the extra push energy is 25.5 MeV. Launching 250000 trajectories with an initial center- 
of-mass energy of 1 MeV above the barrier, we obtained a probability of fusion Pf^g = 
(6.970 ± 0.268) x 10^^^. This was computed using our importance sampling method, with 
an additional force corresponding to A = 1.9/i; about 88% of the trajectories evolving 
under the modified Langevin equation went to fusion. Using Eq.^, our result gives an 
efficiency gain of = 3.5 x 10^! We cannot compare our estimate of Pfus directly to an 
estimate obtained from simulating with the original Langevin equation, since we would need 
to run ~ 10^^ trajectories to have a decent chance of observing even a single fusion event. 
Importance sampling is indispensible in this case: we could not have calculated Pfus using 
direct simulations. 



IV. AN EXACTLY SOLVABLE MODEL 

In addition to the numerical results of the previous section, it is instructive to consider 
an exactly solvable model. Consider a particle in two dimensions which falls at a uniform 
rate Vz from a height h, while experiencing random "kicks" in the horizontal direction: 

z = -Vz , x = iit) , (47) 

where ^(t) represents white noise corresponding to a diffusion constant of unit magnitude: 
{^{t)^{t + s)) = S{s). We assume that the initial horizontal location is zero, i.e. a;(0) = 0, and 
are interested in the horizontal location of the particle when it hits the "ground" {z = 0), 



" This was obtained by running a very large number of simulations for different values of energy 
above the barrier, and then fitting the results to an exponential multiplied by a second-order 
polynomial. The expected error associated with the curve itself is everywhere smaller than the 
smallest of the error bars shown in Fig.^. 
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at time r = h/vz- The motion in the x-direction is a Wiener process, whose solution is a 
Gaussian distribution with a variance growing hnearly with time. Let us suppose we are 
interested in the probabihty that the particle will end to the right of the fixed point Xq, i.e. 
x{t) > xq. Since the ensemble distribution of x(r)-values is a Gaussian (with variance equal 
to r), this probability is given in terms of an error function: 



P[x{t) > xo] = F{xo) = ^[l - erf(xo/v^) 



(48) 



For large values of xq, this probability dies off very rapidly, P ~ e~^o/^'^ / xq, therefore very 
many simulations would be needed to compute P to some desired relative accuracy r. 

Let us now consider an additional force in the form of a constant horizontal "wind", 
pushing the particle in the direction of the point Xq. The modified equations of motion are 
then 



X 



w 



(49) 



where w is the strength of the wind. From Eq.^, it follows that A A = wx(r) — w^r/2. From 
Eq.|27|, we can then compute the minimum number of simulations necessary to obtain P to 
a relative accuracy r, using importance sampling, for a particular value of wind strength: 



N* = — 



1 f „2^F(xo + wr) 



F2(xo) 



1 



(50) 



The efficiency gain is then: 



rpG — ' u 



N* 



F{xo)-F\xo) 



e"^F(xo + wr) - F'^{xo) 



(51) 



For the case Xq = 3.0, r = 1.0 (P = 1.35 x 10~^), we have plotted efficiency gain as a 
function of wind strength, in Fig.^. We see that the optimal wind strength (at which we 
obtain maximal efficiency) is Wopt = 3.157. For this value of w, the number of trajectories 
needed to estimate P using direct simulation of the original Wiener process, is about 220 
times the number needed to estimate P (to the same degree of accuracy) with importance 
sampling. 

It is interesting to note that even for xq = (for which half the trajectories fall to the 
right of Xq), we gain efficiency by using importance sampling. In this case, for r = 1, the 
maximal efficiency gain (about 1.8) is achieved with a wind strength w = 0.6125, as shown 
in Fig.^ 

The maximal efficiency gain grows as the probability P becomes small. Using an asymp- 
totic expression for the error function, we have 



F{xo) 



from which it follows that 




exp(— XQ/2r) (r fixed , xq oo). 



N* 



Xo 
^2 



(xo — wr)^/2r 



(52) 



(53) 
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This formula nicely encapsulates the dramatic efficiency gain achieved for small values of P 
(i.e. large Xq). Without importance sampling, i.e. setting w= 0, the number of trajectories 
needed grows exponentially in x^: Nq ~ exp(xo/2r) (dominant contribution). However, 
using importance sampling, with the optimal wind value (y/opt = Xq/t), the number needed 
grows hnearly with xq: N* ~ xq. Thus, 

K ~ (54) 

Even when Nq is tremendously large, may remain modest (for w = Wopt)- For instance, 
for Xq = 6 and r = 1 we have P[x{t) > xq] = .99 x 10^^. One would need ~ 10^^ direct 
simulations to compute this probability to 10% accuracy. Using importance sampling, one 
can achieve the same accuracy with fewer than 1000 trajectories. Numerical results (with 
w = 6.0) gave us P[x{t) > 6.0] = (1.04 ± 0.10) x 10^^ after only 700 trajectories. 



V. CONCLUSIONS 

In this paper we have developed a method for computing the probabilities of rare events - 
exemplified by the fusion of two nuclei - for processes described by Langevin dynamics. The 
method, based on the idea of importance sampling, is straightforward to implement, quite 
general, and can lead to a very large increase in computational efficiency. For these reasons 
we believe it represents a very practical tool for using numerical simulations to compute 
small probabilities. Indeed, with our method, we were easily able to estimate a fusion 



probability, within a schematic model of nuclear collisions (see the end of Section |T|), that 
would have been essentially impossible to estimate from direct simulations of the process in 
question. We see every reason to expect similar results when combining the method with 
more realistic semiclassical models of nuclear dynamics. 
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FIGURES 



FIG. 1. Thirty trajectories simulated using the schematic model of nuclear collisions, Eq.45. 
The system is ^'^'^Zr+^'^'^'^Zr, at 0.8 MeV above the interaction barrier. Two trajectories lead to 
fusion; the rest to reseparation. 



FIG. 2. The circle indicates the saddle point region of the potential energy, and the arrows 
show the direction of the extra force chosen to push the system toward fusion. Also shown are two 
(deterministic) trajectories, one ending in fusion, the other in reseparation. 

FIG. 3. Convergence of the estimator Pf^^ with number of trajectories simulated (N), for 
both direct simulation (solid line), and using importance sampling (dashed line). 

FIG. 4. Excitation function computed using direct simulation, with 1000 trajectories for each 
point. (The solid line is an analytical estimate extracted from a much larger number of simulations.) 



FIG. 5. Same as Fig.^, but computed using importance sampling instead of direct simulation. 

FIG. 6. Efficiency gain as a function of wind strength, Eq.^, for xq = 3 {P = 1.35 x 10^'^). 
At the optimal wind value, the gain is around 220. 

FIG. 7. Same as Fig.^ but for xq = 0. Here, P = 0.5, so there would be no problem in 
estimating this probability from direct simulations. Nevertheless, there is an efficiency gain of 
nearly two, when using importance sampling. 
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